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ABSTRACT 

We present a numerical study of rapid, so called type III migration for Jupiter- 
sized planets embedded in a protoplanetary disc. We limit ourselves to the case of 
inward migration, and study in detail its evolution and physics, concentrating on the 
structure of the corotation and circumplanetary regions, and processes for stopping 
migration. We also consider the dependence of the migration behaviour on several 
key parameters. We perform this study using the results of global, two-dimensional 
hydrodynamical simulations with adaptive mesh refinement. The initial conditions 
are chosen to satisfy the condition for rapid inward migration. We find that type III 
migration can be divided into two regimes, fast and slow. The structure of the co- 
orbital region, mass accumulation rate, and migration behaviour differ between these 
two regimes. All our simulations show a transition from the fast to the slow regime, 
ending type III migration well before reaching the star. The stopping radius is found 
to be larger for more massive planets and less massive discs. A sharp density drop 
is also found to be an efficient stopping mechanism. In the fast migration limit the 
migration rate and induced eccentricity are lower for less massive discs, but almost do 
not depend on planet mass. Eccentricity is damped on the migration time scale. 
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1 INTRODUCTION 

The theory of disk-planet interaction aims at finding a phys- 
ical description of the gravitational interaction between gas- 
dust protoplanetary disks and protoplanets. A major ad- 
vance, in comparison with the theories of solar system for- 
mation from 20 years ago, was the realization that proto- 
planets can globally migrate and consequently that plan- 
ets can be found on orbits different from the formation 
sites of planetary cores (or gaseous protoplanets, in an al- 
ternative scenario). T h is conclusio n was first reached by 
iGoldreich fc Tremainel l| 19791 . Il980h and worked out in a 
numb er of investigati o ns su mmarised by iLin fc Papaloizoul 
(|l993h and ILin et all (|2000t ). Theories of type I migration 
(planet embedded in gas) and type II (in a disk gap opened 
by the planet's gravity) for a standard stationary accretion 
disk both predict an inward direction of migration, towards 
the star (Ward 1997). Type II migration, which carries a 
gap-opening planet with the viscous flow of material in the 
disk, is capable of moving it outward provided the disk has 
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loca l density gradients supporting outward viscous spread- 
ing |Lin fc Papaloizou|[l993l ). 

The discovery of ex t ra-solar planets llMavor fc Quelozl 
1 19951 . iMarcv et all I2000L IVogt et all |2002| ) triggered a re- 
newed interest in the migration process. For the so-called 
hot Jupiters, giant planets residing at a distance of typi- 
cally just 0.05-0.1 AU from the host star, an in situ forma- 
tion appears impossible in either the core accretion or the 
disk instability scenario of planet formation. Therefore, in 
the current consensus, hot Jupiters are understood as bodies 
that started forming at orbital radii about 100 times larger, 
where availability of ice enhances the rate of growth of plane- 
tary cores, and later migrated toward the star by interacting 
wit h the disk (for an exte nsive recent review of migration, 
see iPapaloizou et alj|2007l ). Another important issue in the 
theory of planet formation is to understand the intermediate 
(e = 0.2 — 0.5) orbital eccentricities of planets. While in this 
paper we concentrate primarily on migration, our numer- 
ical investigations also produce results for the eccentricity 
evolution of a planet migrating in a disk. 

Until recently, disk-planet interaction was understood 
entirely in terms of just one type of coupling mecha- 
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nism, involving resonances. Disks are perturbed by a planet 
most strongly at Lindblad resonances, which are the gas- 
dynamical equivalents of mean motion resonances in celestial 
mechanics. Disks respond by launching density and bending 
waves. These waves take the angular momentum and energy 
flux deposited by the perturber and transport it away from 
their launching point, to later deposit these quantities in 
the disk. The non-axisymmetric wave pattern acts back the 



planet, providing the t orque and energy flu: 



orbital evolution (e.g., iGoldreich fe Tremain 



to drive its 



19791 . Il980t ). 



The action of Lindblad torques is not only responsible, for 
orbital migration of types I and II, but also for the shep- 
herding effect in the disk-satellite interaction as well as the 
gap opening phenomenon. 

However, in addition to the Lindblad resonances, effects 
often described as corotational resonances can play a signif- 
icant and even dominant r ole, as found independe ntly by 
iMasset fc Papaloizoul |2003l) and lArtvmowiczl (120041 ). Those 
effects are not truly resonant, in the sense of depending on 
strict mean-motion commensurabilities, beyond a long-term 
1:1 commensurability of the planet with gas orbits librat- 
ing with it. We will therefore refer to them as corotational 
torques or corotational flows. The corotational region cor- 
responds to the horseshoe and tadpole orbit region in a 
frame corotating with the planet. It occupies an annulus 
a few Roche lobe or Hill sphere radii away from a circular 
planetary orbit. Gas in this region alternates between being 
closer and further from the star than the planet. In doing 
so, the gas may flow at small density but significant velocity 
through the gap region, if a gap exists. 

Misunderstood for a long time as unimportant (partly 
because of an assumed smooth density distribution in a 
disk, and partly due to previous estimations done for a 
fixed-orbit, i.e., non-migrating, planet, which automati- 
cally eliminates the effect), the effects of the flow of gas 
across the corotational region can provide a far larger 
torque on the planet than the competing, wave- launching 
Lindblad interactions. This leads to a new type of very 
rapid migration (called type III), which apart from speed, 
also distinguishes itself by having a far smaller predilec- 
tion for inward migration than the standa r d type I or 
type II varieties jMasset fc Papaloizoul | 2003l : lArtvmowiczl 
l2004l ; |Papaloizoull2005l ; lArtvmowiczll2006l). This type of mi- 
ratio n was studied numerically by Masset fc Papaloizoul 
2003), who performed two-dimensional simulations of a 
freely moving planet in steady-state migration for a whole 
range of imposed fixed migration rates. Global, high res- 
olution two- and three-dimensional simulations of a freel y 
migrating planet were performed bv lD'Angelo et ail (|2005h . 
iPapaloizoul (120051') considered local shearing box simulations 
in 2-D, and lArtvmowicz fc de Val Borrol <|2007f > in 3-D. 

This paper is the second one in a series on numer- 
ical simulations of pl anet migration. In the first paper 
ijPeplinski et al. I l2007al . henceforth Paper I) we showed how 
two new, physically motivated, ingredients in our model, 
namely the use of a modified local-isothermal approxima- 



1 We shall neglect the mention of energy flow and only describe 
the torque or angular momentum flow. Only the latter drives 
migration of bodies whose eccentricity does not grow (a typical 
result of numerical simulations such as ours). 



tion and a correction for self-gravity effects, are necessary 
to achieve consistent and convergent numerical results. In 
the current paper, using these methods, we study the evo- 
lution of type III migration and the physical dependencies 
between the model parameters and the outcome (mostly, mi- 
gration rate and range, and the flow of gas around and onto 
the planet). In other words, we focus on the physics rather 
than numerics of type III migration. We limit ourselves in 
this paper to models describing the inward migration of a 
freely evolving planet. Paper III (IPeplinski et al.ll2007bl ) will 
describe the fast outward migration, which is not a simple 
mirror case of the inward migration. 

The layout of the paper is as follows. In Section [5] we 
outline the mechanism of type III migration. Sections [3] 
and [3] briefly summarise the basic model equations and the 
numerical method. In Section [S] we describe and analyse in 
detail our 'standard case', treating its orbital evolution, evo- 
lution of torques, and the flow pattern near the planet. Sec- 
tion [6] deals with the dependence of the migration process 
on various parameters. In Sections [7] and [8] we analyse the 
stopping of type III migration and the eccentricity evolution. 
Finally, in Section [9] we summarise our findings. 



2 TYPE III MIGRATION REGIME 

As outlined above, type III migration is mostly driven by 
the flow of material through the co-orbital region. Although 
the corotation torque dominates, the differential Lindblad 
torque still retains some importance by modifying the space 
of allowed solutions for migration rates and generally bias- 
ing the migration to become inwardly directed, an in type I 
migration. The source of the corotation torque is the inter- 
action between the planet and orbit-crossing fluid elements. 
The fluid exchanges angular momentum with the planet 
when it executes a U-turn at the end of each horseshoe- 
like streamline. In the absence of planet migration, dissi- 
pation or a time-dependent effect such as an initial tran- 
sient, the fluid elements in the co-orbital region circulate 
on closed orbits (in the corotating frame), resulting in a 
zero net corotation torque in a smooth disc. This cancella- 
tion is often called the saturation of the corotational torque. 
In this stationary case only an actively maintained den- 
sity or more generally, a specific vorticity gradient in the 
planet's vicinity can produce a non-zero corotation torque. 
The density variation can be due to the internal disc struc- 
ture (e.g., a n edge of a dead zo ne, where viscosity changes 
abruptl y, cf. lArtvmowiczl l2006lh or to disc heating by the 
planet dPaardekooper fc Mellemj|2 006) or due to MHD tur- 
bulence (|Nelson fc Papaloizoull2004) . 

For a migrating planet, the situation is clearly more 
susceptible to a non-cancellation of torques. The horseshoe 
region becomes asymmetric and some streamlines are no 
longer closed (those belonging to the disk). On the other 
hand, a region of librating, closed orbits forms. It is carried 
along together with the migrating planet as 'ballast' if it is 
denser than the disk, or as a 'balloon' if it is less dense. In 
the latter case it so to say, propels the planet along its orbit 
by the imbalance of gravitational pull of gas in front and in 
the back of the planet. In fact, the initial direction of type III 
migration being toward the denser side of the disk (to com- 
pensate for the flow of mass and angular momentum of gas 
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from that dense to the rarefied side, and conserve the angu- 
lar momentum of planet-disk system globally) , migration al- 
ways results in the librating orbit region being under-dense, 
i.e. the second case just mentioned. We can thus describe 
the disk with a planet in type III migration as having at 
least a partial gap, but one which is front-back asymmetric. 

This c onfiguration provides a positi ve feedback for the 
migration. iMasset fc Papaloizoul (|2003T > showed that the 
corotation torque is proportional to the migration rate d 
for slow migration (defined quantitatively below). At large 
imposed migration rates, they found that the torque dimin- 
ishes slightly after reaching a maximum, and obtained the 
migration speed estimate of 



2IY 



«a(M£ - 5m) ' 



(1) 



where Flr is the differential Lindblad torque, Mp — Mp + 
Mr is all of the mass content within the Roche lobe (sum 
of the planet and gas mass), and 5m is the co-orbital mass 
deficit: the difference between the mass that the horseshoe 
region would have if it had a uniform surface density equal 
to the upstream surface density E s and the actual horseshoe 
region mass Mhs, 



5m = 47raa; s E s — Mhs = 47raa; s (E s — E g ), 



(2) 



where E g is an average density in the horseshoe region (gap 
region jfj. Equation ([1} clearly cannot be applied to the most 
interesting case, where the denominator of Eq. JTJ tends to 
zero, the mass deficit becomes comparable with or larger 
tha n M^, and we expect a tr a nsition to fast migration. 

lArtvmowicz fc Peplinskil (|2007T) considered a simple 
model in which the angular momentum of a fluid element 
undergoes a sudden change at orbital conjunction with the 
planet. They found a critical migration speed a = df at 
which one of the ends of the horseshoe-like orbit cannot 
reach conjunction and dis-attaches from the planet, becom- 
ing a tadpole orbit. This speed is the ratio of the half- width 
of the corotational horseshoe region x a (the subscript 's' 
stands for separatrix distance) and one-half of the libration 
time (synodic period) on the critical horseshoe/tadpole or- 
bit. The df corresponds to the maximum corotation torque 
in this model and provides a convenient fiducial speed unit 
for type III migration. For a Keplerian disc with angular 
speed at planet's semi-major axis O, df has the form 



at = 



3x 2 s fl 
8ira ' 



(3) 



In terms of this fiducial speed, we call migration with 
|d|/df < 1 slow migration, and migration with |d|/df > 1 
fast migration. 

To simplify the notation, we define two useful, non- 
dimensional quantities. The first is the speed of migration Z 
that we will extensively use to describe a momentary state 
of our numerical models, 



z = e. 

df 



(4) 



Z expresses the ratio of the migration timescale Tmigr (de- 
fined as x B /a) to the libration timescale Tub- We therefore 
distinguish between the fast (\Z\ > 1) and slow migration 
limits (\Z\ < 1). 

The second dimensionless quantity is Ma, expressing 
the co-orbital mass deficit 



Mx 



5m/ 



47rax s (E s — E g 
Ml 



(5) 



lArtvmowicz fc Peplinskil (|2007h derived an analytical 
theory for the whole range of migration rates based on a 
guiding centre approximation. They obtained x s = 2.47 » 
2.5 Hill sphere radii (Roche lobes) of the planet, and con- 
cluded that this value provides a good approximation to the 
numerically determined estimators of the corotational re- 
gion's half-width. The corotation torque is a linear function 
of Z for Z < 1 (slow migration limit) and becomes constant 
for Z > 1 (fast migration limit): 



5m 



I^r = sign(d)— df (l - [max(0, 1 - \Z\ 



0,a. 



(6) 



Unlike IMasset fc Papaloizoul (|2003h . 

lArtvmowicz fc Peplinskil (|2007h do not predict migra- 
tion to have the runaway or self-accelerating character of 
an instability. Rather, a steady migratio n speed Z results 
from the torque balance equation (cf. also IPapaloizou et all 
120071 ) in which the following three torques are in balance: 
the torque required for migration with speed Z according 
to orbital mechanics of nearly circular orbits, the corotation 
torque Fcr, and the Lindblad resonant torque Flr. This 
allows one to define Z as a function of Ma and determine 
the space of allowed solutions for a given value of Tlr. The 
relation between Z and Ma is linear in the fast migration 
limit \Z\ — 2Ma/3 (for Ma > 1.5) and becomes nonlinear 
in the slow migration limit (Ma < 1.5). 

From Eqs. © and we can see that For ~ Mp in 
the fast migration limit and increases slower than linearly 
in the slow migration limit. This means that the planet's 
orbital evolution should be weakly dependent on the planet 
mass in the fast migration limit. However, in the slow migra- 
tion limit any increase of the planet mass should decrease 
the migration rate significantly. We verify these predictions 
below. 

To conclude this description of type III migration we 
want to point out that it in general depends on the migration 
history, because the way the horseshoe region is populated 
depends on the previous evolution. In the case of strong 
variations of d the streamlines of the horseshoe region are 
not exactly closed, so the co-orbital mass deficit can be los t 
and the fast migration would stop l|Papaloizou et al.ll2007l ). 



2 A note of caution: despite its name and physical units, 5m pro- 
vides a measure of surface density deficit rather than the actual 
mass difference of the librating region against the background 
disk, owing to the variable shape of that region (5m becomes the 
actual mass deficit only if migration is infinitesimally slow and 
the librating region is an annulus of a constant width 2x B ). 



3 DESCRIPTION OF THE PHYSICAL MODEL 

The full description of the disc model and the numeri- 
cal method was given in Paper I. In this section we give 
only a short description of the governing equations and our 
methodology. 
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3.1 Disc model 

We adopt in our simulations a two-dimensional, infinitesi- 
mally thin disc model and use vertically averaged quanti- 
ties. We work in the inertial reference frame, in the Carte- 
sian coordinate system (x,y,z). The plane of the disc and 
the star-planet system coincides with the z = plane. The 
centre of mass of the star-planet system is initially placed 
at the origin of the coordinate system. 

The gas in the disc is taken to be inviscid and non-self- 
gravitating. The evolution of the disc is given by the two- 
dimensional (x,y) continuity equation for E and the Euler 
equations for the velocity components v = (v x ,v y ) 



0, 



9Ev 



•V(Euv) +VP = -EV$, 



(7) 
(8) 



where P is two-dimensional (vertically integrated) pressure, 
and $ is the gravitational potential generated by the proto- 
star and planet 



$ = $s + $p = 



GM S 
r - rs 



GAR 



rp 



(9) 



The positions and masses of the star and the planet we de- 
note by rs, Ms, rp and Mp respectively. Mp is the effective 
planet mass, as explained below. The gravitational poten- 
tial close to the star and the planet is softened according to 
Eq. 5 from Paper I, which contains a parameter r so ft, the 
smoothing length. 

We do not consider the energy equation, but use a local 
isothermal approximation which prescribed the temperature 
as a function of both distance to the star and to the planet 
such that the global disc has a scale height h s with respect 
to the star, and the circumplanetary disc a scale height h p 
with respect to the planet. This equation of state was called 
EOS2 in Paper I. 

In Paper I we also showed that it is necessary to include 
some sort of self-gravity in the neighbourhood of the planet. 
For this we use the extra acceleration given by Eq. 13 from 
Paper I for the gas near the planet. The size of the region 
to which this correction is applied is characterised by a pa- 
rameter r cnv . For this correction we introduced the concept 
of the effective planet mass, Mp , which we normally take to 
be Mp — Mp + Af so ft , the sum of the planet mass and the 
gas mass within the softening radius. For experiments where 
we want to explicitly ignore the mass accumulation near the 
planet, we take Mp = Mp. 



3.2 Equation of motion for the star and the planet 

The goal of the current study is to investigate the orbital 
evolution of the planetary system due to the gravitational 
action of the disc material. Since the calculations are done 
in the inertial reference frame, the equation of motion for 
the star and the planet have the simple form: 



rs = - 



rp = 



GMp(r s -r P ) 
\rp - r s \ 

GMsjrp - r s ) 
rp - rsl 3 



G(r s - r)dM D (r) 

i |3 ' 

\r - r s \ 

G(r P - r)dM D {r) 
\r — rpl 3 



(10) 



(11) 



In both cases the integration is carried out over the disc 
mass Mb included inside the radius r^ac- 



3.3 Setup description 

In the simulations we adopt non-dimensional units, where 
the sum of star and planet mass Ms + Mp represents the 
unit of mass. The time unit and the length unit are chosen 
to make the gravitational constant G = 1. This makes the 
orbital period of Keplerian rotation at a radius a — 1 around 
a unit mass body equal to 2n. However, when it is necessary 
to convert quantities into physical units, we use a Solar-mass 
protostar Ms = Mq, a Jupiter-mass protoplanet Mp = 
Mn,, and a length unit of 5.2AU. This makes the time unit 
equal to 11.8/27T years. 

In all simulations the grid extends from —4.0 to 4.0 in 
both directions around the star and planet mass centre. This 
corresponds to a disc region with a physical size of 20.8 AU. 

3.3.1 Initial conditions 

The initial surface density Emit is given by a modified power 
law: 



E ta 



V>(r c )E (r c /r ) c 



(12) 



where r c = |r — rc| is the distance to the mass centre of the 
planet-star system, vq is a unit distance, and if> is a function 
that allows introducing a sharp edges in the disc (see Fig. 2 
in Paper I). 

The disc mass is described by the disc to the primary 
mass ratio 



E in it(ro)7rro S 7rro 



(13) 



M s M s 

In the simulations fip> ranges from 0.0025 to 0.005. For the 
Minimum Mass Solar Nebula (MMSN) /id = 0.00144 (for 
qe = —3/2). We investigate different density profiles by 
changing as from —1.5 to —0.5. 

To enforce rapid migration the planet is introduced in- 
stantaneously. Since we focus on type III migration only and 
do not analyse the problem of orbital stability inside a gap, 
we do not introduce the planet smoothly nor keep it on a 
constant orbit for the time needed to create a gap. For the 
simulated cases of inward migration the initial density gra- 
dient given is sufficient to start rapid migration. 

In most simulations the initial planet mass is such that 
Mp/Ms = 0.001 (i.e. one Jupiter mass, M^, for a one-solar- 
mass star). We investigate the migration of the planets with 
Mp/Ms equal 0.0007 and 0.0013 too. The planet always 
starts on a circular orbit of semi-major axis equal 3.0. 

The aspect ratio for the disc with respect to the star is 
fixed at h s = 0.05, whereas the circumplanetary disc aspect 
ratio hp lies in the range 0.4 to 0.6. 

The smoothing length of the stellar potential, r so f t , is 
taken to be 0.5. For the planet this parameter was cho- 
sen be a r so f t = 0.33-Rh, where the Hill radius is given by 
Rh = a[M P /(3M s )] (1/3) . The characteristic size of the re- 
gion where the self-gravity correction is applied is r cnv = 
0.5.Rh. 

In our simulations we use an outflow-inflow boundary 
condition with a so-called killing wave zone next to the 
boundaries, where the solution of the Euler equations are 



© 0000 RAS, MNRAS 000,1111221 



Type III migration; Inward migration 5 



smoothly connected with the initial disc in sub-Keplerian 
rotation. For a full description see Paper I, Sect. 2.3.2. 

To track the flow of material through the corotation 
region, we use a tracer fluid that has a value of 1 in the 
(initial) corotation region (atnjt — 2Rh) < r < (dinit + 2Rh), 
and zero outside it. This allows us to distinguish between 
gas captured by the planet in the horseshoe region and gas 
flowing through the corotation region. 



4 NUMERICAL METHOD 

For our simulations we use the Flash hydrodynamics code 
version 2.3 written by the FLASH Code Group from the ASC 
/ Alliance Center for Astrophysical Thermonuclear Flashes 
at the University of Chicagqj in 1997. Flash is a modular, 
adaptive-mesh, parallelised simulation code capable of han- 
dling general compressible flow problems. It is designed to 
allow users to configure initial and boundary conditions, 
change algorithms, an d add new physics mo dules. It uses 
the Paramesh library l|MacNeice et al.l l2000h to manage a 
block-structured adaptive mesh, placing resolution elements 
only where they are needed most. 

For our purpose, the code is used in the pure hydrody- 
namical mode in two dimensions, and the refinement crite- 
ria are set to achieve high resolution around the planet. Our 
simulations use a lowest resolution mesh of 800 cells in each 
direction, and a square region around the planet is refined. 
The maximal cell size in the disc (lowest level of refinement) 
is 0.01, the minimal cell size is 0.00125 (4 levels of refine- 
ment, in practice corresponding to 1.8% of the minimum Hill 
sphere radius). 

Further details about the numerical method are given 
in Paper I. 



5 INWARD MIGRATION - STANDARD CASE 

In this section we describe the inward migration of our stan- 
dard case: a Jupiter-mass planet Mp = 0.001 with the cir- 
cumplanetary disc aspect ratio h p — 0.4 in a disc with 
Ho = 0.005 and qe = —1.0; the disc has no outer edge. The 
surface density is taken to be initially constant inside the 
gravitational softening for the star, giving a small jump in 
Ei n it at r = 0.5. However, this jump is too small and located 
too close the star to influence the planet's migration. In the 
simulation the effective gravitational mass of the planet was 
increased by the mass content within the smoothing length. 

5.1 Orbital evolution 

The planet's orbital evolution is shown in Fig. [J|. Looking 
at the semi-major axis (upper left panel) we can divide the 
planet's orbital evolution into two stages. The first one is the 
rapid migration stage, which lasts for about 60 orbits and 
corresponds to the type III migration regime. During this 
stage, the planet migrates rapidly due to the co-orbital flow 
and is not able to clear a full gap. The migration rate a (mid- 
dle left panel) is initially close to constant around —0.0065, 

3 http://fiash.uchicago.edu 

4 All plots in this paper present data averaged over 5 periods. 



but after t — 40 orbits, |d| starts to decrease rapidly. This 
behaviour is caused by the strong mass accumulation in the 
planet's proximity (see Sect. 15.31 and Fig. |3). As the migra- 
tion slows down some of this mass is lost again, giving a close 
to zero near t = 55 orbits. The dimensionless migration rate 
Z (middle right panel) shows essentially the same behaviour 
but since at ~ a~" ,5 (Mp/Ms) 2 ' 3 , it decreases even when a 
is almost constant. The important transition from the fast 
{\Z\ > 1) to the slow (\Z\ < 1) migration limit takes place 
at about 47 orbits. 

The migration slows down as a = 1 is approached. 
In the second stage (after 60 orbits) the planet migrates 
slowly enough to allow a gap to open up. To avoid confu- 
sion with the concept of the slow migration limit (defined as 
\Z\ < 1), we will refer to this second stage as the gap open- 
ing stage. The |d| decreases slowly and is equal 1.5 x 10~ 4 
at t = 100 orbits, giving a migration time-scale tm = a/d 
equal 1.3 x 10 4 yr. Note that this is not classical type II mi- 
gration, as the total torque is not dominated by the the dif- 
ferential Lindblad torque. However, as the corotation torque 
|F| slowly decreases with time, the system may ultimately 
make the transition to type II migration. 

The time evolution of the total torque F exerted by the 
gas on the planet is plotted in the lower left panel. In the 
rapid migration stage the migration is driven by the corota- 
tion torque Fcr and the total torque F is dominated by Tcr, 
the differential Lindblad only playing a negligible role. Un- 
like the migration rate, |T| grows during the first 40 orbits, 
and then rapidly drops due to the strong mass accumulation 
in the planet's proximity and the transition from the fast to 
the slow migration limit. The corotation torque is closely re- 
lated to the migration rate, as illustrated in the lower right 
panel. The plot shows V versus the non-dimensional migra- 
tion rate Z. What is most notable is that the total torque 
|T| grows linearly with \Z\ for the slow migration limit. In 
the fast migration limit the relation between Z and T is non- 
linear and the torque saturates around Z ~ —1.5. |F| reaches 
it s maximum value at Z ~ —1 .7. It is consistent with Fig. 9 
in lMasset fc Papaloizoul (120031 ). 

In the gap opening stage T evolves toward zero (it is 
equal -1.8 x 10~ 4 and -1.1 x 10" 4 at 75 and 100 orbits 
respectively), but is still larger than the predicted differen- 
tial Lindblad torque. This is caused by a small mass outflow 
from the interior of the Roche lobe, which increases |r|. In 
models where mass accretion on the planet is allowed, this 
mass outflow is much more limited and the final torque is 
an order of magnitude smaller. 

The upper right panel in Fig.[T]shows the changes of the 
eccentricity e. The eccentricity grows rapidly during the first 
few orbits and reaches 0.03, after that decreasing with time. 
In the gap opening stage e becomes e < 0.01. We postpone 
further discussion of the eccentricity evolution to Sect. [8] 



5.2 Flow structure in the co-orbital region 

In this section we discuss the relation between the coro- 
tational torque Tcr and the asymmetry of the horseshoe 
region. The latter is connected to the the non-dimensional 
migration rate Z, since Z expresses the ratio of the migra- 
tion to libration time scales. In the slow migration limit 
(\Z\ < 1) T m i gr is longer than Tub and the horseshoe region 
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Figure 1. Orbital evolution for the standard case. The upper left and upper right panels show the evolution of the planet's semi-major 
axis a and eccentricity e. The evolution of the migration rate a and the non-dimensional migration rate Z are presented in middle left 
and middle right panels. The lower row shows the total torque T exerted by the gas on the planet. Left and right panels present the 
torque as a function of time and as a function of the non-dimensional migration rate Z, respectively. 



can extend the full 2-7T in azimuth. For \Z\ > 1 the horse- 
shoe region shrinks into a single tadpole-like region. Once 
this happens the rest of the co-orbital region is taken up by 
a flow from the inner to the outer disk. We will denote this 
flow as well as the region it goes through as the co-orbital 
flow. 

In the simulation Z evolves from —3 to (middle right 
panel in Fig. [TJ , and thus the morphology of the co-orbital 
region should show an evolution from tadpole to a regular 
horseshoe shape. To study this we use the tracer fluid to 



follow the gas that was placed in the initial corotational 
region (ainit — 2Rh < r B < ainit + 2Rb). Figure [2] shows 
the evolution of this tracer fluid, with the top left panel 
displaying the initial condition. 

The top right and middle left panels present the mass 
fraction during the fast migration phase for Z = —2.35 
(t = 20 orbits) and Z = -1.4 (t = 40 orbits). Due to the 
planet's fast inward migration, most of the gas placed ini- 
tially in the co-orbital region is left in its initial position and 
only a small fraction of mass is captured by the planet in 
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Figure 2. The mass fraction of the gas that initially was placed in the corotational region fljnjt — 2/2jt <C fs <C c^init "t - where Ts is 

the distance to the star and a; n ; t is the planet's initial semi-major axis. The mass fraction goes from (gas from the inner disc only) up 
to 1 (gas from corotation only). The top left panel shows the initial condition. The top right (t = 20 orbits, Z PS —2.35) and a middle 
left (t = 40 orbits, Z pa —1.4) present the fast migration limit. The middle right (t = 50 orbits, Z —0.4) and a bottom left (t = 53 
orbits, Z ft: —0.05) show the slow migration limit. The last panel displays the mass fraction during the gap opening stage (t = 80 orbits, 
Z « -0.046). All plots are in a co-moving reference frame with the planet at (x,y) = (3,0), (2.25,0), (1.44,0), (1.1,0), (1.05,0) and 
(0.99, 0) respectively. 
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the horseshoe region. The gas flow there is strongly asym- 
metric and the region shrinks to a single tadpole-like region 
around the L4 libration point. It has relatively sharp edges 
and the gas captured by planet almost does not mix with 
the gas crossing the co-orbital region. The azimuthal extent 
of the horseshoe region grows slowly with decreasing \Z\ and 
is about 7r and 1.2tt for Z equal —2.35 and —1.4 respectively. 
This means that the co-orbital flow takes up a wide range in 
azimuth, that is only weakly dependent on Z. Consequently, 
the mass flux of the gas crossing corotation depends mostly 
on the initial density profile, and the co-orbital torque Tor 
is a linear function of Ma , almost insensitive to the changes 
of Z. We should keep in mind that Ma is not a really an 
independent variable, since it depends on the way the horse- 
shoe region is populated, which depends on the migration 
history and thus indirectly on the previous values of Z. 

The volume of the co-orbital region is given by Vcr = 
4nax a , and thus decreases for inward migration, as can be 
seen in Fig. [2] This is important since the mass deficit is 
a function of Vcr and of the difference in the density flux 
between the co-orbital flow (E s ) and the horseshoe region 
that moves together with the planet (E g , also see Eq.[2]). For 
a constant value of E s — E g , the torque would thus decrease 
as Vcr, i.e. as a 2 . A slow down of inward migration can thus 
be caused by this simple geometrical reason. We discuss this 
and other stopping mechanisms in Sect. [7] 

During the entire phase of fast migration, there is a 
small amount of the gas leaving the horseshoe region. This 
is visible in Fig. [2] as a light-blue leading spiral. 

As \Z\ drops, the horseshoe region grows in azimuth and 
the co-orbital flow region shrinks. For \Z\ ~ 1 the horseshoe 
region fills the whole co-orbital region, and the simulation 
enters the slow migration limit. In the presented simula- 
tion the surface density at 5.2 AU is equal 523 g cm -2 and 
E ~ r" 1 , and the transition between fast and slow migration 
limit happens when the planet is at 6.2 AU. In the slow mi- 
gration limit phase the flow through the corotation is limited 
to a narrow stream at the boundary of the Roche lobe and 
the asymmetry in the horseshoe region is relatively small. 
However, this small asymmetry is sensitive to small changes 
in Z, giving a strong dependence of Tor on Z. As the flow 
becomes more symmetric both globally (shape of the horse- 
shoe orbits) and locally (shape of orbits in the planets vicin- 
ity; see Sect. 15. 3| ), T diminishes rapidly. 

The middle right and bottom left panels present the 
mass fraction for the slow migration limit for Z = —0.4 
{t — 50 orbits) and Z — —0.05 (t = 53 orbits). In this phase 
the planet's radial motion is slow and a gap is being cleared. 
Although the azimuthal extent of the horseshoe region is 
close to 27r for both values of Z, the shape of the horseshoe 
streamlines differ between the two snapshots. In the earlier 
case an asymmetry is still visible, but in the later one some 
gas from the inner disc has been trapped in the tadpole 
region around the L5 Lagrange point and the asymmetry is 
starting to disappear. In this phase of migration the second 
tadpole region is created and the gas from the co-orbital flow 
mixes with the gas carried along with the planet. On both 
panels there is a visible flow through the co-orbital region, 
but in the second case this flow is dominated by gas leaving 
the planet's vicinity. This lowers V and gives a « at about 
55 orbits. 

The last panel displays the mass fraction after the 



planet has almost opened a gap (t = 80 orbits, Z — —0.046). 
In this case the asymmetry of horseshoe streamlines is al- 
most invisible, and the gas flow through the co-orbital region 
is very small. We can still see a difference between the tad- 
pole regions around the L4 and L5 Lagrange points. Most of 
the gas around L4 has been carried along with the planet, 
but the gas around L5 was mostly captured during the last 
phase of migration. 



5.3 Mass accumulation and flow structure in the 
planet's vicinity 

In the previous section we described the relation between the 
non-dimensional migration rate and the asymmetry of the 
gas flow in the full co-orbital region. We will now focus on 
the local gas flow near to the planet. Just as we saw for the 
whole co-orbital region, the planet's fast radial motion also 
strongly influences the gas evolution inside its Roche lobe, 
causing the circumplanetary disc to acquire an asymmetric 
shape. The co-orbital flow enters the Roche lobe and after 
interacting with the gas in the planet's proximity, will partly 
be captured into the circumplanetary disc. We have found 
this mass accumulation to be an important ingredient in 
our simulations of type III migration, determining the flow 
structure in the planet's vicinity and the torque generated 
within the Roche lobe. 

In our standard simulation, the effective gravitational 
planet mass was increased by the mass content within the 
smoothing length Mp = Mp + M so f t , where the smooth- 
ing length correspond s to the "surface" of the planet 
l|Peplinski et al.l l2007al ). The time evolution of Mp is pre- 
sented in the left panel of Fig. [3] During the first 40 or- 
bits the migration rate is approximately constant, as is the 
mass accumulation rate. During this time the mass con- 
tained within the smoothing length M ao f t reaches about 45% 
of the initial planet mass. The mass of the gas inside the 
whole Hill sphere amounts to about 55% of the initial mass, 
so most of this material is found inside the planet's gravita- 
tional smoothing length. 

The mass accumulation rate increases sharply when the 
planet starts to slow down at t = 40 orbits. In the relatively 
short time of 10 orbits M so f t increases up to 90% of the 
initial planet's mass and a strong pressure gradient forms 
around the planet. This pressure gradient is supported by 
the gas inflow into the circumplanetary disc. The large mass 
accumulation rate is made possible by our use of a locally 
isothermal equation of state that makes the circumplanetary 
disc aspect ratio independent of the gas density. For the 
same reason the accumulated mass leaves the Hill sphere 
shortly after the rapid migration has stopped. This affects 
the total torque T giving the positive value of a during this 
short period (Fig. [T] middle panels) . 

When the simulation reaches the second, gap opening 
stage, M ao ft remains approximately constant at 56% of Mp. 

The relation between the planet's effective mass Mp 

and Z is presented in right panel of Fig. [3] Mp is low and 
almost constant during the fast migration limit phase and 
rapidly grows when Z reaches —1.7, causing the total torque 
to decrease (lower right panel in Fig. [TJ. The next impor- 
tant change takes place in the slow migration limit phase 
at Z ~ —0.6, when the planet starts to open a gap and 
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Figure 3. The evolution of the effective planet mass for the standard case. The left and right panels present the effective planet mass 
as a function of time and of the non-dimensional migration rate, respectively. The initial planet mass is Mp = 0.001. During the first 

40 orbits the migration rate d is almost constant and the rate of the mass accumulation Mp remains constant too. Mp increases when 
migration slows down. After the planet enters the gap opening stage at about 60 orbits the planet loses about 14% of its mass and at 
the end of simulation Mp = 0.00156. 



Mp starts to decrease. We note that the evolution of Mp is 
strongly dependent on choice for the circumplanetary disc 
aspect ratio h p , but for the given value of h p — 0.4 the evo- 
lution of Mp is almost independent of the other parameters, 
such as r so f t . For more discussion on this, see Sect. [5] and 
Paper I. 

Since the mass accumulation rate determines the struc- 
ture inside the Roche lobe, we now discuss the flow in the 
planet's vicinity for four different stages of the simulation: 
the rapid migration stage with and without the fast mass 
accumulation in the circumplanetary disc, the stage of mass 
outflow from the planet's proximity and the final gap open- 
ing stage with an approximately constant effective planet 
mass. They are presented in Fig. [4] The plots show the 
surface density and the flow lines. The flow lines close to 
the border of different regions have been made more visi- 
ble. There are five regions marked in the plots: inner (I) and 
outer (II) disc, the co-orbital region (III and IV) and the 
circumplanetary disc (V). It should be realized that the gas 
evolution inside the Roche lobe is variable and we present 
just a few snapshots to describe important differences be- 
tween different stages of simulation. 

The upper left panel of Fig. [4] corresponds to the fast 
migration limit with Z ~ —2.35 and a relatively low mass 

accumulation rate Mp. In this case region III is the horse- 
shoe region and IV the co-orbital flow. The orbits in the 
circumplanetary disc are found to be strongly asymmetric, 
compressed at the side of the co-orbital flow and stretched 
at the side of the horseshoe region. The asymmetry of the 
bow shocks and the flow lines in the co-orbital region is 
also visible. The whole horseshoe region is shifted backwards 
with respect to the planet's radial motion, and the left bow 
shock is compressed into a straight line, whereas the right 
bow shock is stretched and more curved. In this phase the 
mass accumulation is still rather modest; closer investigation 
shows the circumplanetary disc to be fed by both regions 
III and IV. The averaged density and the torque asymme- 
try distribution for this stage of migration are presented in 
Fig. El 



The stage of strong mass accumulation near the transi- 
tion from the fast to the slow migration limit (Z ~ —1-1) is 
presented in upper right panel in Fig. [4] The asymmetry of 
the bow shocks and the flow lines in the co-orbital region is 
still visible, however the orbits in the circumplanetary disc 
have become less regular. The interior of the Roche lobe is 
strongly disturbed by the gas inflow from the inner disc I 
and the horseshoe region III and the flow lines are very time 
variable. In spite of this, the circumplanetary disc no longer 
shows the asymmetry seen in the previous phase, which ex- 
plains the lower value for the torque F. 

The asymmetry of the bow shocks and the circumplan- 
etary disc disappears in the next stage. The lower left panel 
shows the flow lines near the planet during the slow mi- 
gration limit phase (Z « —0.05) and strong mass outflow. 
As we saw in the previous section, in this stage the horse- 
shoe region extends over the whole azimuthal range and the 
planet starts to open a gap. The whole co-orbital region 
also becomes more symmetric and regular. The gas leaves 
the circumplanetary disc flowing along the bow shocks and 
influences the flow lines at the border of the Roche lobe, but 
the orbits inside the circumplanetary disc are regular. 

The last panel (lower right) presents the stage of the gap 
opening planet \Z\ <C 1, with an almost constant mass in the 
circumplanetary disc. In this stage the gap is almost cleared 
and the co-orbital region is symmetric. Although the mass 
inside the Hill sphere is almost constant some outflow from 
the Roche lobe can still be seen. The most important differ- 
ence with the previous stages is the shape of the flow lines 
in the planet's vicinity. In the previous stages the gas was 
orbiting the planet creating a circumplanetary disc, however 
in this last stage the large pressure gradient prevents circular 
orbits in the planet's proximity. We have to note that this 
stage looks different for models with an accreting planet (see 
Paper I, Sect. 5.1.2). In that case the planet is surrounded 
by the Keplerian-like sub-disk and there is no gas outflow 
from the Roche lobe during the last phase of migration. Thi s 
is consistent with what was found bv lD'Angelo et all (2003). 

After having described the flow structure we now con- 
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Figure 4. The surface density and flow lines in the planet's vicinity for the standard case at four different stages of migration. The 
upper left panel corresponds to the fast migration limit with Z fs —2.35, t = 20 orbits. The upper right and lower left panels show 
slower migration {Z m —1.1, t = 46 orbits and Z «s —0.05, t = 56 orbits) with the circumplanetary disc quickly gaining and losing mass 
respectively. The last panel (lower right) presents the stage of the gap opening planet (\Z\ <1, t= 100 orbits) with an almost constant 
mass of the circumplanetary disc. The colour scale is logarithmic. The plotted domain is square region of the size of 6i?u- The flow lines 
close to the border of different regions are made more visible. There are five regions marked in the plots: inner (I) and outer (II) disc, 
the co-orbital region (III and IV) and the circumplanetary disc (V). On the upper panels the regions III and IV denote the horseshoe 
and co-orbital flow regions, respectively. 



sider the torque distribution in the planet's vicinity. We fo- 
cus on the first stage, i.e. rapid migration with a low mass ac- 
cumulation rate (i = 20 orbits) . This corresponds to the up- 
per left panel of Fig. [4] In the left panel of Fig. [5]we present 
the surface density averaged over 5 orbits for Z ~ —2. The 
plot covers 5Rn and the planet is placed at the centre. The 
red line shows the Roche lobe size. The asymmetry between 
the left and right bow shock which we described above, is 
clearly visible. Although the left bow shock is sharper and 
has a higher density close to the planet, the overall amount 
of mass at the right bow shock is higher, since the density 
increase close to the right bow shock extends over a much 
wider area. The plot also shows the averaged surface density 
to be lower in the horseshoe region (positive values of y) than 
in the co-orbital flow (negative values of y). This asymme- 
try causes the right bow shock to have a larger contribution 
to the total torque T than the left one. This is further il- 



lustrated in the right panel of Fig. [S] where we show the 
absolute value of the sum F c (x p ,y p ) + r c (:Ep, — y p ), where 
r c is the torque generated by a single grid cell and (x p ,y p ) 
is the position of the cell with respect to the planet. The 
borders of the regions giving positive (+) and negative (-) 
contributions to the total torque are shown with pink lines. 
The left and the right bow shocks give positive and nega- 
tive contributions, respectively. Although the absolute value 
of the torque sum is similar for both shocks, the right bow 
shock is much wider and its contribution to F dominates. 
The region between the two shocks gives a negative contri- 
bution showing that the torque generated by the co-orbital 
flow is stronger than the torque generated by the horseshoe 
region. 

The interior of the Roche lobe shows three extremes. 
The two outer ones give a positive contribution and the mid- 
dle one a negative. These extrema are the result of the asym- 
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Figure 5. The density distribution and the asymmetry in the torque distribution in the planet's vicinity for the fast migration limit 
with the low mass accumulation rate (Z fa —2.35, t = 20 orbits). The left panel shows the surface density averaged over 5 orbits. The 
right panel presents the absolute value of the sum r c (x p ,y p ) + r c (iEp, — y p ). The borders of the regions giving positive (+) and negative 
(-) contributions to the total torque are shown by thick pink lines. The red line shows the Roche lobe size. The plots cover 5-Rh and the 
planet is placed at the centre. The colour scale is logarithmic. 



metry of the circumplanetary disc, which as we saw above, is 
compressed at the side of the co-orbital flow and stretched 
at the side of the horseshoe region. However the sum of 
these extremes nearly cancels, and the overall torque gen- 
erated within the Roche lobe is found to be negative. This 
once more illustrates the point made in Paper I, that the 
interior of the Roche lobe cannot be treated as a system iso- 
lated from the co-orbital flow and that the torques from this 
inner region have to be taken into account when studying 
type II I migration (for more discussion, see IPeplinski et ahl 
l|2007al )). 

We would like to repeat here that we only show one 
snapshot here, but that the gas flow in the planet's vicinity 
is highly variable, with changes in the shapes of bow shocks 
and the circumplanetary disc, causing the large variations 
in r seen in the lower panels of Fig. [T] 
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5.4 Surface density 

To gain further inside into the migration process, we now 
study the surface density distribution in the disc. A colour 
scale plot of the global surface density distribution for the 
fast migration limit phase at t — 20 orbits (Z fa —2.35) is 
displayed in Fig. [5] In this phase the planet does not yet open 
a gap due to its fast radial motion. For the same reason the 
surface density in the inner disc does not show significant 
differences from the initial condition. There are two lower 
density regions at the planetary radius corresponding to the 
horseshoe region (positive y) and the co-orbital flow region 
(negative y). The detailed structure of the surface density 
distribution in the planet vicinity is shown in Fig. [5] 

To study the time evolution of the surface density E we 
use a series of plots showing this quantity in the radial and 
azimuthal direction. Fig. [7] presents the azimuthal average 
of E over 2-7T at times t — 0, 10, 30, 50 and 60 orbits (curves 
1, 2, 3, 4 and 5 respectively). The average surface density is 



Figure 6. Global surface density for the standard case at t = 20 
orbits. The planet has migrated from a = 3 to a = 2. The colour 
scale is logarithmic. The planet's location is (x,y) = (2.3,0) 



weakly modified outside the planet's co-orbital region. Dur- 
ing the rapid migration phase (curves 1, 2, 3 and 4) no gap 
is cleared, but the planet is followed by a lower density re- 
gion. The main reason for this depletion in E in the planet's 
neighbourhood is the mass accumulation in the planet vicin- 
ity and the fact, that after an interaction with the planet the 
gas crossing the co-orbital region is spread over a region with 
a radial size of a few x s . The last curve in Fig. [7] corresponds 
to the gap opening stage (d approaching the value for type 
II migration). 

Fig. [8] shows E along an azimuthal cut through the 
planet position for time t = 10, 30, 50, 60 and 100 orbits 
(curves 1, 2, 3, 4 and 5 respectively). It illustrates the dif- 
ference between the density of the horseshoe region (E g ) 
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Figure 7. The azimuthal average of the surface density £ for 
times 4 = 0, 10, 30, 50 and 60 orbits (curves 1, 2, 3, 4 and 5 
respectively). During the rapid migration phase no gap is formed 
(curves 1 to 4). The last curve corresponds to the gap opening 
stage The planet position is visible as a spike in the density profile. 




together with the azimuthal average of the surface density 
(curve 2). Curve 1 shows the initial density profile. The up- 
per left and right panels show these quantities during the 
fast migration limit phase (t = 30 and 40 orbits, respec- 
tively) and the lower left and right panels the same during 
the slow migration limit phase (t = 50 and 60 orbits). The 
cut through L4 gives the approximate width of the co-orbital 
region and the density E g inside the horseshoe region. Simi- 
larly the cut through L5 gives the approximate value of the 
density of the gas crossing the co-orbital region E s . The dif- 
ference between E a and E g is the main reason for the rapid 
migration. 

In the fast migration limit the planet carries along the 
low density horseshoe region which has well defined, sharp 
edges and a width x w m 2Rn- E g is lower than the initial 
density at the initial planet position r = 3 and almost does 
not change in time. In contrast, E s is close to the initial 
density profile at the current planet's position. E g and E s 
differ from each other in the whole co-orbital region. The 
cut through L5 shows also the low density region following 
the planet. It is caused by the spreading of the gas over a 
few x w and compression at the two shocks in the outer disc 
region (see Fig. [6} . 

In the slow migration limit (lower left panel) the den- 
sity around L5 drops significantly, as the horseshoe region 
gets more extended and the co-orbital region becomes more 
symmetric. At the same time the density inside the horse- 
shoe region increases slowly. The asymmetry between the 
density cuts through L4 and L5 is now only visible at the 
inner edge of the co-orbital region. In the last, gap opening 
stage of migration Z <C 1 (lower right panel) the co-orbital 
region becomes symmetric and the density cuts through L4 
and L5 are similar. The planet starts to opens a gap and the 
density in the horseshoe region decreases. 



Figure 8. The surface density £ along an azimuthal cut through 
the planet position for times t = 10, 30, 50, 60 and 100 orbits 
(curves 1, 2, 3, 4 and 5 respectively). The negative and positive 
values of <j> correspond to the co-orbital flow and horseshoe region, 
respectively. The planet is located at <j> = 0. 



and the co-orbital flow region (Es), that drives type III mi- 
gration. The region of co-orbital flow and the horseshoe re- 
gion are found at the negative and positive values of (f>, re- 
spectively. The fast migration limit phase is represented by 
curves 1 and 2. Initially the asymmetry in E is only located 
inside the Roche sphere (only barely visible on the plot due 
to adopted scale), but in time it grows in azimuth show- 
ing a large difference in E between the L4 and L5 Lagrange 
points. During the fast migration limit phase the L5 point is 
moved in the direction of the planet (compare the position 
of the density maximum for positive (j> in curves 2 and 5). 
The asymmetry decreases in the slow migration limit phase 
(curve 3) and becomes unimportant when a gap starts open- 
ing up (curves 4 and 5). During the whole simulation the 
surface density has maxima at the L4 and L5 points and 
decreases in the direction of the planet. 

To further illustrate the asymmetry, Fig.[9]shows radial 
surface density cuts through L5 (curve 3) and L4 (curve 4), 



6 DEPENDENCE ON THE SIMULATION 
PARAMETERS 

After having described type III migration in detail for 
one standard case, we now consider the dependence of the 
planet's orbital evolution on the most important simulation 
parameters. 



6.1 Planet's mass 

The first parameter we consider is the planet's mass. In our 
investigation we concentrate on the orbital evolution of giant 
planets. The results of the simulations with the initial mass 
Mp equal 0.0007 and 0.0013 (curves 1 and 3 respectively) 
are presented in Fig. 1101 together with our standard case of 
M P = 0.001 

The upper left and upper right panels show the evolu- 
tion of the planet's semi-major axis a and migration rate a. 
During the first 40 orbits all systems evolve almost exactly 
the same way, even though the effective planet masses (lower 
right panel) differ significantly. This period corresponds to 
the fast migration limit \Z\ > 1 for all three models. In this 
phase the orbital evolution is independent of the planet mass 
since Z ~ Ma, and both the non-dimensional migration rate 



Z and the mass deficit Ma are proportional to M p 
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Figure 9. The evolution of the density profiles in the standard case. Each panel shows the initial density profile (curve 1), the azimuthal 
average of the surface density (curve 2; the planet position is visible as a spike) and the surface density cuts through the Lagrange points 
L5 (curve 3) an L4 (curve 4). The panels upper left, upper right, lower left and lower right show the density profiles after 30, 40, 50 and 
60 orbits, respectively. 



Once the planets enter the slow migration limit, the re- 
lation between Z and Ma becomes more complicated and 
their evolution starts to differ. In this phase, |d| decreases 
rapidly and the planet is gradually settles into a type II 
like migration. Since Z diminishes faster for more massive 
planets, this transition to the slow migration limit happens 
sooner for more massive planets, and as a result more mas- 
sive planets stop their type III migration at larger orbital 
radii. The final planet's semi-major axes are 0.75, 0.95 and 
1.11 for the initial mass M P equal 0.0007, 0.001 and 0.0013, 
respectively. 

The relation between the planet mass and its orbital 
evolution is further illustrated in the lower left panel, where 
the non-dimensional migration rate Z is plotted as a func- 
tion of the planet's position in the disc. A larger planet mass 
corresponds to a smaller \Z\ and thus displaces the orbital 
evolution curve to the right. 

The lower right panel shows the effective mass evolu- 
tion. For a constant circumplanetary disc aspect ratio h p , 
higher mass planets accumulate a larger amount of gas in 
the circumplanetary disc. The mass accumulation rate is 
relatively low during the fast migration limit and increases 

with decreasing \Z\. There is a visible increase of Jlfp for 
Z ~ —1.7 in all three cases. Similarly the mass outflow 
from the Roche lobe starts at Z m —0.6 independently of 



the planet mass (and thus happens earlier for more mas- 
sive planets). The final planet masses are 0.00095, 0.00155, 
0.0023 for M P equal 0.0007, 0.001 and 0.0013 respectively. 
This corresponds to a relative increase of 35%, 55% and 76% 
after 100 orbits respectively, and thus the relation between 
Mp and Mp is not linear. 

The lowest mass planet shows some evolutionary fea- 
tures not seen in the other two. First, its mass accumulation 
rate shows two instances of rapid growth, the first one at 
about 27 orbits (Z w —2.5), and the second one at 45 or- 
bits (Z ~ —1.7). Second, after 80 orbits the lowest mass 
planet shows oscillations in its migration rate. After 120 or- 
bits these oscillations have damped out, and the migration 
rate d becomes similar to that of the more massive planets. 
These oscillations are caused by a change of shape of the 
flow lines in the planet's vicinity from regular orbits in the 
circumplanetary disc to a non regular flow modified by the 
strong pressure gradient. 

6.2 Effects of circumplanetary disc aspect ratio 

The next parameter we consider is the circumplanetary disc 
aspect ratio h p . As shown in Paper I, this parameter corre- 
sponds to the temperature of the circumplanetary disc, and 
determines how much material can be accumulated there. 
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Figure 10. Results of the simulations for different planet masses. Curves 
0.0007, 0.001 (standard case) and 0.0013, respectively. The upper left and 
major axis a and the migration rate d^The non-dimensional migration rate 
evolution of the effective planet mass Mp are presented in the lower left and 
limit \Z\ > 1) all systems evolve almost exactly the same way even though 
to differ after reaching the slow migration limit. Even though the migration 
migration limit, the non-dimensional migration rates Z differ slightly, since 
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Time [Orbits] 

1, 2 and 3 correspond to the initial planet mass Mp equal 
upper right panels show the evolution of the planet's semi- 
Z as a function of the planet's position in the disc and the 
lower right panels. During the first 40 orbits (fast migration 
the (effective) planet mass may vary. The simulations start 
rate a is almost^identical for all simulations during the fast 
df depends on Mp. 



Since the mass accumulation influences type III migration, 
it is interesting to study the effect of this parameter. 

In addition to our standard case of h p — 0.4, we ran 
simulations using values of 0.5 and 0.6. We do not use a 
value lower than 0.4, since the tests presented in Paper I 
showed that a model with 0.3 does not fully converge for 
the resolution used in our simulation:^. 

The results for all three cases (labelled 1, 2 and 3) 
are presented in Fig. [11] We use the same set of plots as 
in Fig. 1101 the evolution of the planet's semi-major axis a 
and the migration rate d in the upper panels, and the non- 
dimensional migration rate Z as a function of the planet's 
position in the disc together with the evolution of the effec- 
tive planet mass in the lower panels. 

As expected, changing h p has a strong effect on the evo- 
lution of the effective planet mass Mp. As we saw above, the 
standard case (0.4, curve 1) shows a phase of strong mass ac- 
cumulation, followed by mass loss from the circumplanetary 
environment after migration has slowed down. The other 

5 The results for h p = 0.3 are actually a more extreme version of 
the 0.4 results, see Paper I. 



two cases do not show such strong mass accumulation, and 
the higher h p , the less mass is accumulated. In fact curves 
2 and 3 show that mass outflow already starts at t ~ 20 
orbits, i.e. during the rapid migration stage. 

Considering the evolution of a and a we notice some 
small differences. For h p = 0.5 and 0.6 the flow asymmetries 
in the planet's vicinity during the fast migration limit phase 
become less pronounced, leading to a lower total torque 
This results in a somewhat lower value for the migration rate 
a during this phase. The start of mass outflow at t ~ 20 or- 
bits in these two models leads to a small increase of the mi- 
gration rate. Consistent with the results from the previous 
section, the lower final values of Mp cause the two thicker 
disc cases to migrate closer to the star. Their type III migra- 
tion stops at a ~ 0.7. It is noticeable that the evolution of 
the system in the slow migration limit is almost independent 
of the circumplanetary disc aspect ratio for h p ^ 0.5. 

Tracing the evolution in the Z — r plane (lower left 
panel), all three simulations show similar behaviour until 
Z ~ —2 is reached (lower left panel). After \Z\ drops below 
1.7 the model with h p — 0.4 starts rapid mass accumulation 
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Figure 11. Results of the simulations for the different circumplanetary disc aspect ratios h p . Curves 1, 2 and 3 correspond to h p equal 
0.4, 0.5 and 0.6 respectively. The panels are the same as in Fig. 1101 



in the circumplanetary disc, causing it to evolve differently 
from the other two models. 

We have to keep in mind that the assumption of a con- 
stant circumplanetary disc aspect ratio is a rather crude one, 
and in reality this parameter will react to the mass accretion 
into this disc. More realistic calculations would either need 
to treat the thermal evolution of the circumplanetary disc 
self-consistently, or use some approximate recipe relating h p 
to Z and the mass accumulation rate. However, the results 
in this section show that large changes in h p only have a 
limited effect on the migration behaviour of the planet. 

6.3 Total disc mass 

In the previous two sections we discussed parameters con- 
nected to the planet and its circumplanetary disc. In this 
section and the next we consider parameters defining the 
circumstellar disc, namely the total disc mass (given by the 
disc to the primary mass ratio /xd) and the slope of the 
initial density profile (described by the exponent as), see 
Sect. I3XT1 

To study the relation between the total disc mass and 
the planet migration we compare a simulation with /id = 
0.0025 to our standard case of 0.005. The other parameters 
are identical between these two simulations. The results are 
presented in Fig. [12] 

As can be expected the planet migrates slower and less 



far in the lower mass disc. The stage of the rapid migra- 
tion lasts longer (85 orbits instead of 50) and the transition 
between the fast and slow migration limits happens at 60 or- 
bits, instead of 40. After the rapid migration ends the planet 
is locked in the disc at a = 1.3 (a — 0.95 for the standard 
case). 

There is no simple relation between the total disc mass 
and the position where the planet is locked into the disc, or 
the duration of the rapid migration stage. This is because 
the migration rate in type III migration depends both on 
the co-orbital mass deficit Ma and on the asymmetry in 
the co-orbital region. Changing the total disc mass changes 
the mass inside the co-orbital region, but the asymmetry 
depends on the current and previous values of the migra- 
tion rate. However, /xd is an important parameter, since the 
amount of the mass in the co-orbital region defines the re- 
gion where rapid migration can take place. Since the volume 
of the co-orbital region is proportional to a 2 , for the stan- 
dard smooth density profile with as ^ —1.5 it is easier to 
start r apid migration at larger orb ital radii (see the discus- 
sion in lMasset fc Papaloizoul [20031. The radii r = 0.95 and 
r = 1.3 are the outer borders of the region, where rapid mi- 
gration for a Jupiter-mass planet cannot be started for /xd 
equal 0.005 and 0.0025 respectively. 

The initial value of the migration rate is given by the 
amount of the gas in the co-orbital region. After about 10 or- 
bits the non-dimensional migration rates are Z ~ —2.7 and 
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Figure 12. Results of the simulations for different total disc masses. Curves 1 and 2 correspond to the disc to the primary mass ratio 
/iQ equal 0.005 (standard case) and 0.0025 respectively. The upper left and right panels show the evolution of the semi-major axis and 
the effective mass; the non-dimensional migration rate Z as a function of time and of radius are presented in the lower left and lower 
right panels. 



Z ~ — 1.4 (lower left panel in Fig. [12}, their ratio indeed 
matching the ratio of the disk masses. Later both systems 
evolve differently. The non-dimensional migration rate de- 
creases quickly during the whole stage of rapid migration 
for the more massive disc, but for the lower mass disc \Z\ 
increases during the first 35 orbits and then decreases much 
slower. More similarities between both simulations are visi- 
ble after plotting Z as a function of the planet's position in 
the disc (lower right panel in Fig. I12|l . Neglecting the first 
25 orbits (region outside r « 1.9 and r « 2.5 for the high 
and low mass disc respectively) the shape of both curves is 
similar, but the curve for the more massive disc is shifted 
inward by 0.3 in radius. This number corresponds exactly to 
the shift of the boundary of the region where rapid migration 
is allowed. Decreasing the disc mass thus has a similar effect 
as increasing the planet mass (lower right panel in Fig. I10[) . 

The effective planet mass (upper right panel in Fig. I12p 
depends on the total mass of the disc and on the non- 
dimensional migration rate. It grows slowly in the fast mi- 
gration limit phase until \Z\ drops below 1.7 after which the 
mass accumulation rate increases rapidly for both models. 
The maximum value of the planet mass for both models is 
reached at Z ~ —0.6 at which point the mass outflow from 
the circumplanetary disc starts. The values are 1.56Mp and 
1.22M P . 



6.4 Initial density profile 

The next parameter is the exponent as- The standard case 
has as = — 1 (curve 2 in Fig. I13p . and we performed two 
additional simulations using —0.5 and —1.5 (curves 1 and 3). 
All other parameters were taken identical to the standard 
case. Because of our definition of the surface density, the 
disc densities are equal at r — 1, but differ at the initial 
planet position r = 3. The as = —0.5 case has the largest 
initial mass in the co-orbital region, and therefore evolves 
the fastest. It reaches a « 1 after about 43 orbits. For the 
as = —1.0 and —1.5 cases the planet settles at the same 
position after about 51 and 64 orbits respectively. Because 
this position is near r — 1, all three models agree on the 
final orbital radius, since this is where the surface density is 
identical between the models, and at a value that is too low 
to support further rapid migration. 

The exponent as determines the average rate with 
which the non-dimensional migration rate is decreasing and 
influences the relation between Z and the planet's position in 
the disc, which can be seen in the lower right panel. The evo- 
lution of the systems differs during the first 16 orbits only. 
Later on the relation between Z and a becomes independent 
of as ■ However, as we will see below this is mostly caused by 
the evolution of the effective planet mass which is presented 
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in the upper right panel. In the fast migration limit Mp de- 
pends on the initial mass of the co-orbital region (largest for 
«s = —0.5) and grows slowly. The mass accumulation rate 

Mp increases when the planet makes the transition to the 
slow migration limit. The increase is steep for as = —0.5 
and relatively smooth for as = —1.5. What is important 

is that Mp has similar values in all presented models and 
the highest value of Mp is reached in the simulation with 
steepest density profile, since the mass accumulation phase 
takes the longest time there. The final mass of the planet is 
equal 1.62M P , 1.56M P and 1.72M P for as equal -0.5, -1.0 
and —1.5 respectively. 

To test how the mass accumulation near the planet in- 
fluences the results we ran a set of simulations where the 
correction for self-gravity is turned off by keeping the planet 
mass constant at Mp = Mp = 0.001. In these simulations we 
also introduced an inner disc edge (see Paper I Fig. 2). This 
edge is not significant for the subject discussed in this sec- 
tion, but it is relevant for the discussion on stopping mech- 
anisms given in the next section. The inner gap starts at 
r = 1.3 and the density is constant at E = 0.00016 inside 
r = 0.85. 

The results of these simulations are shown in Fig. 1141 
The upper left panel shows the planet's orbital evolution. 
In all simulations the planet stops its rapid migration at 
the disc edge. Similarly the transition between the fast and 



the slow migration limit takes place at the same position, 
a — 1.28 close to the initial density maximum. The interac- 
tion with the disc edge also leads to a temporary outward 
directed migration, although this effect was found to be sen- 
sitive to the resolution used. The stopping of rapid migration 
at the disc edge is caused by the loss of flow asymmetry in 
the co-orbital region due to the interaction with the low den- 
sity region. A more detailed discussion is presented in the 
next section. 

Like in the previous model the initial amount of the gas 
within the co-orbital region determines the planet's orbital 
evolution and the disc with the smaller density gradient gives 
faster migration. In a similar way it influences the initial 
mass of the gas inside the Hill sphere. The later evolution 
of Mh depends on as- For as = —0.5 Mh remains almost 
constant until the planet reaches the disc edge, and then 
starts to decrease. For the other models this mass remains 
approximately constant until Z > —1.7, and then starts to 
increase. When the planet reaches the disc edge Mh has a 
similar value for all models, and then decreases with the 
same rate. 

The relation between the non-dimensional migration 
rate and as is illustrated on lower right panel showing Z 
as a function of the planet's position in the disc. Here we 
see that without mass accumulation, the curves do not over- 
lap. In the fast migration limit Z is approximately constant 
for as = —1.5 and changes fastest for as = —0.5. This could 
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Figure 14. Results of the simulations for the different initial disc profiles with constant Mp = Mp = 0.001. Curves 1, 2 and 3 correspond 
to as equal —0.5, —1 and —1.5 respectively. The layout of the panels is the same as in Fig. 1121 except that here the upper right panel 
shows the evolution of the mass inside the Hill sphere. 



be due to the fact that Z is related to the specific vorticity 
w — (A x «)/£, which has an almost constant value in the 
disc for as = —1.5 and decreases with the radius for bigger 
as . Alternatively, this relation may also be explained by the 
fact, that in the fast migration limit Z ~ Ma ~ a 2 (E s — E g ). 

All the curves merge at Z » —2 and r ~ 1.45, since 
E s at becomes similar in all simulations and the interaction 
with the disc edge starts to be important. Later on all the 
systems evolve in a similar way due to the interaction with 
the disc edge, however the temporary outward migration is 
stronger for the model with as = —0.5. 



7 STOPPING TYPE III MIGRATION 

In this section we address the processes that stop rapid type 
III migration. We consider two possible causes. In the first 
one migration slows down due to the decrease of the volume 
of the co-orbital region, and in the second one the planet 
stops at an inner disc edge. When discussing stopping of 
type III migration we mean a radical slow-down of migration 
to a rate similar to type II migration, or \Z\ <C 1. However, 
an exact stopping, Z = 0, can be observed in some cases. 

Since the non-dimensional migration rate Z is a func- 
tion of the co-orbital mass deficit Ma, to stop a type III 
migration we have to diminish Ma. This can be achieved 



by either decreasing the mass of the co-orbital region or 
by removing the asymmetry in this region. The migration 
rate depends on the amount of angular momentum trans- 
ferred between the planet and the gas crossing the planet's 
orbit. To start and maintain rapid migration the impulse 
that the planet gets from the orbit crossing gas has to be 
strong enough to move the planet into the undisturbed re- 
gion of the disc where it can interact with a new portion of 
gas. This impulse depends on the change in the radial po- 
sition of the gas element crossing the orbit given by :r s and 
the asymmetry in the co-orbital region given by E s — E g 
(Eq. [2J . Assuming maximal asymmetry and putting E g = 
leaves us with two parameters, E s and a; s , whose product is 
proportional to the mass of the co-orbital region. This pa- 
rameter allows us to qualitatively understand the stopping 
mechanism for type III migration and the dependence of Z 
on the disc total mass and on as. It is however difficult to 
get quantitative results, since E g in practice is not zero, and 
depends on the migration history. 

Our simulations show that fast migration is not allowed 
in the region where the mass of the co-orbital region Mcr 
is smaller than the planet mass0 Since for as > — 2 this 



6 The exact calculations show, that at the stopping radius Mcr 
ranges between 1.8 and 2.3 of Mp taken at the stopping time. 
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Figure 15. The evolution of the density profiles during the planet migration in a disc with an inner edge. Each panel shows the initial 
density profile (curve 1), the azimuthal average of the surface density (curve 2; the planet position is visible as a spike) and the surface 
density cuts through the Lagrange points L5 (curve 3) an L4 (curve 4). The panels upper left, upper right, lower left and lower right 
show the density profiles after 40, 50, 60 and 70 orbits, respectively. 



parameter decreases with decreasing distance to the star, 
it is easier to start type III migration further out in the 
disc. The migration rate decreases with distance to the star 
and thus there will be an inner region where fast migra- 
tion cannot occur. This means that in most cases the planet 
will not be able to migrate inward all the way by type III 
migration. The value of the radius where the planet stops 
rapid migration depends on the surface density profile and 
the planet mass. This radius decreases with increasing 
and decreasing as- Planets with smaller masses should be 
able to migrate further in, but with a lower migration rate 
since df ~ Mp 2 ^ 3 . Interestingly, this analysis implies that 
the situation is reversed in discs with steep density profiles, 
qe < —2. In this case it is easier to start fast migration close 
to the star and the inward migrating planet should increase 
the migration rate with decreasing distance to the star. 

Mcr decreases mostly due to the shrinking of the ra- 
dial extent of the co-orbital region. This can be seen in 
Fig. [9] where the the azimuthal average of the surface den- 
sity (curve 2, and the surface density cuts through the L5 
(curve 4) an L4 (curve 5) points are plotted. The planet 
moves through the disc leaving the surface density profile 
slightly changed (curve 1 is the original density profile). Af- 
ter 30 and 40 orbits the strong asymmetry between L5 and 
L4 that is driving migration is visible. However, this asym- 



metry disappears after 50 orbits even though the mass dif- 
ference between the inner and outer edge of the co-orbital 
region still exists. The migration stops since the jump in the 
radial position of gas crossing the co-orbital region is too 
small to remove a substantial amount of angular momen- 
tum from the planet. 

The other mechanism for stopping type III migration is 
the interaction with a steep density change, such as a disc 
edge. At a significant density jump Mcr will diminish con- 
siderably, and asymmetry in the co-orbital region will disap- 
pear. Figure [15] shows the evolution of the density profiles in 
the simulation presented at the end of Sect. 16.41 The upper 
left, upper right, lower left and lower right panels show the 
density profiles at t = 40, 50, 60 and 70 orbits respectively. 
At these times Z » -2, -1.5, -0.05, 0.0. At t = 70 orbits 
the planet is locked in the disc and gap opening has started. 
The upper left and upper right panels show the density dis- 
tribution during the fast migration limit Z < — 1 and are 
similar to the upper left and upper right panels in Fig. [9] 
The lower left panel shows the planet interacting with the 
disc edge. At this time the mass difference between the in- 
ner and outer edge of the co-orbital region disappears and 
the planet gets a positive value of Z. However we found this 
temporary outward migration to be sensitive to the grid res- 
olution (decreasing for higher resolution). 
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The results show that a Jupiter is massive enough to de- 
stroy the initial density jump and fill the inner empty region. 
That is why when the planet finally settles, the inner and 
outer disc have almost equal surface densities. This actually 
prevents the planet from stopping completely. The situation 
is different for lower mass planets that cannot change signif- 
icantly the density profile. They can be trapped close to the 
position of maximum density for a long time and grow slowly 
l|Masset et al.ll2006i ). To have this happen to a Jupiter-mass 
planet would require an additional mechanism that removes 
material from the inner disc and keeps the radial density gra- 
dient sharp. Such a mechanism would keep the planet at a 
constant orbit due to the gas flow from the co-orbital region 
until the horseshoe region has completely cleared. However, 
if the rate of mass accretion onto the star is too small to fill 
the co-orbital region and keep a constant gas flow through 
the co-orbital region the planet would still move away from 
its position at the disc edge. 

In paper III we will encounter two additional mecha- 
nisms that can stop type III migration, but which are more 
evident in the outward migration case. 



8 ECCENTRICITY EVOLUTION 

In this section we discuss the eccentricity evolution and its 
dependency on different simulation parameters. The results 
of different models are presented in Fig. 1161 see also the top 
right panel of Fig. Q]for our standard case. 

The left panel of Fig. [TS] shows the time evolution of 
the eccentricity e for the model with the constant effective 
planet's mass Mp = Mp — 0.001 (curve 1) and the model 
including the gas accretion on the planet (curve 2). In the 
second model the planet's mass was increased by the gas 
mass removed from the planet's proximity (for more discus- 
sion see Paper I Sect. 5.1.2). It shows that during the rapid 
migration phase (lasting for about 70 orbits) the evolution of 
the eccentricity is almost independent of the planet's mass 
and the mass of the circumplanetary disc. The eccentricity 
grows rapidly during the first few orbits and reaches 0.018, 
and diminishes with time on the timescale of migration or 
faster. The differences between the models show up in the 
slow migration phase, where the first model loses gas from 
the Hill sphere and the eccentricity starts to oscillate due 
to these gas motions in the circumplanetary disc. In the 
second model accretion has removed most of the mass of 
the circumplanetary disc, and the eccentricity slowly decays 
without oscillations. We found similar results in simulations 
with the Mp = Afp+M so f t and different initial planet masses 
Mp. In all of these cases the time-scale for the eccentricity 
damping during the rapid migration phase is independent 
on the planet mass. We found it to be weakly dependent 
on other simulation parameters such as the density profile 
exponent as- 

However, the time-scale for the eccentricity damping 
does depends on the total disc mass. The right panel in 
Fig. [16] shows the time evolution of e for the models with 
/io equal 0.005 and 0.0025 (curves 1 and 2). The value of 
e after a few first orbits depends on j^u and is about 0.03 
and 0.017 for the first and the second model respectively. 
In both models the eccentricity is damped during the rapid 
migration phase, and the corresponding damping time-scale 



is about 33 and 69 orbits. This suggests that the damping 
time-scale is proportional to 1/Vd 

A probable explanation for the eccentricity damping 
is the presence of dense gas in the region of the strong e- 
damping co-orbital Lindbl ad resonance s . Alth ough the ec- 
centricity driving theory of lArtvmowicd l| 19931 ) does not ex- 
plicitly account for the modification of gas flow in the corota- 
tional region and is not meant for gap-opening planets, it can 
provide an approximate upper limit to the strength of eccen- 
tricity damping by density wave emission at the Lindblad 
resonances. Eccentricity driving is a result of the compe- 
tition between external Lindblad resonances and co-orbital 
Lindblad resonances, the latter being about 3 times stronger 
than the former. Since type III migration implies that the 
gap around the planet is asymmetric, let us assume that the 
mean gas density in the co-orbital resonance region near the 
planet is reduced by a factor of 2 with respect to the case of 
a smooth disk, considered by I Artvmowic z ( 1993), equivalent 
to the simplified view of zero gas density in the libration re- 
gion and unperturbed values elsewhere. Eccentricity damp- 
ing occurs in direct proportion to the gas density. Therefore, 
the resulting eccentricity damping timescale should in our 
standard case be larger than Artymowicz's estimate by a 
factor of (3 - l)/(3/2 - 1) = 4. In a disk with fj, D = 0.005, 
h s = 0.05 everywhere, and qe = —3/2, we obtain the fol- 
lowing lower limit on the damping timescale at a = 5 AU: 

e/e > (0.008//i)yr (14) 

(where fi — 0.001 for Jupiter), giving 8 years (about one 
orbit) at Jupiter's radius. This is much shorter than the 
numerically determined timescale for e to drop from e = 
0.02718 to e = 0.01, which equals about 30 periods, as can 
be seen in Fig.[T] We conclude that the analytical description 
is consistent with our calculations but does not provide a 
good guidance. A theory better adjusted to the corotational 
gas flow streamlines and/or eccentric corotational torques is 
needed. It is interesting to note that eccentricity damping 
and migration can operate on a similar, short timescales. 



9 CONCLUSIONS 

In this paper we presented a study of inward directed type 
III migration of Jupiter-mass planets embedded in a disc. 
We used two dimensional numerical simulations, performed 
in a Cartesian coordinate system and the inertial reference 
frame. We focused on the detailed flow structure in the 
planet's vicinity and used an adaptive mesh refinement to 
achieve high resolution inside the Roche lobe. We used a 
modified version of the usual local-isothermal approxima- 
tion, where the temperature depends on the distance to both 
the star and the planet, and also added a correction for the 
gas self- gravity near to the planet. 

Type III migration is driven by the co-orbital torque 
Tcr exerted on the planet by the gas crossing the co-orbital 
region. The value of this torque depends on the mass flux 
crossing the planet's orbit and on the density asymmetry in 
the co-orbital region. This asymmetry can be characterised 
by the non-dimensional migration rate Z = a/ at which gives 
the ratio between the the migration timescale (defined as 
Xs/a) and libration timescale. We can thus divide type III 
migration into a slow (\Z\ < 1) and a fast (\Z\ > 1) migra- 
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Figure 16. The time evolution of the eccentricity e. The left panel shows the results of the models with the constant effective planet's 
mass Mp = Mp = 0.001 (curve 1) and the model including the gas accretion on the planet (curve 2). The right panel presents the results 
of the models with different total disc masses. Curves 1 and 2 correspond to fip, equal 0.005 (standard case) and 0.0025 respectively. 



tion limit. The shape and the structure of the horseshoe re- 
gion differ between these two limits as the azimuthal extent 
of the horseshoe region decreases with \Z\. In the slow mi- 
gration limit the horseshoe region fills the whole co-orbital 
region and the flow asymmetry is limited to the planet's 
direct vicinity. Both tadpole regions are present, although 
their shape changes with Z and the positions of the station- 
ary points move away from the Lagrangian points L4 and 
L5. In the fast migration limit the whole horseshoe region 
shrinks to a single tadpole-like region. 

In the fast migration limit the planet's radial motion is 
too fast to allow a gap to be opened (Fig. [6] and [7]). Instead 
there is a strong gas flow through the corotation region, 
moving gas from the inner disc (the region with r < a — x s ) 
to the outer disc (the region with r > a + x s ). In the inertial 
reference frame the planet together with the gas captured in 
the horseshoe region moves fast with respect to the gas that 
was initially placed in the inner disc (Fig[2} . Note that there 
is no global motion for most of the gas in the disc in this case. 
Although the fluid captured by the planet in the horseshoe 
region moves together with the planet, the gas crossing the 
corotation changes its initial position in the disc only by 
order of x s . This is an important difference with type II 
migration, where the planet is locked in the disc and moves 
on the viscous time scale. In the slow migration limit, when 
the migration timescale becomes longer than the libration 
timescale, the planet starts clearing a gap, and gradually 
makes the transition to type II migration. 

Another difference between the slow and fast migration 
limits is the relation between Fcr and Z. The torque ex- 
erted on the planet depends mostly on the structure of the 
flow lines in the horseshoe region and in the circumstellar 
disc. In the slow migration limit the flow through the coro- 
tation is limited to the a narrow stream at a boundary of 
the Roche lobe and the asymmetry of both the circumplan- 
etary disc and the horseshoe region is relatively small. This 
asymmetry is sensitive to the small changes in Z, giving a 
strong dependence of the total torque T (dominated by the 
corotational torque) on Z. On the other hand in the fast 
migration limit the asymmetry of the circumplanetary disc 
and the horseshoe region is significant, the gas can cross the 



co-orbital region over a wide azimuthal range and the to- 
tal torque is only weakly dependent on Z. This results in 
a linear grow of the total torque |T| with \Z\ in the slow 
migration limit and a saturation of |P| around Z ~ —1.5. In 
the fast migration |F| decreases slowly with \Z\. 

The torque also depends on the rate of mass accumu- 
lation in the planet's vicinity, which is related to the non- 
dimensional migration rate (Fig. [3}, and on the assumed 
scale height for the circumplanetary disc. For h p = 0.4 the 
effective planet mass grows slowly in the fast migration limit, 
since in this case the horseshoe region contracts into a sin- 
gle tadpole-like region, and most of the fluid crossing the 
co-orbital region does not enter the Roche lobe. The gas or- 
bits in the circumplanetary disc are found to be strongly 
asymmetric, and the gas flow into the Roche lobe is limited. 
The asymmetry of the horseshoe region and the circumplan- 
etary disc diminishes with decreasing \Z\, and the mass ac- 
cumulation rapidly increases for Z » —1.7. This causes the 
torque generated inside the Hill sphere to drop significantly 
and shortly after the planet enters the slow migration limit, 
where the co-orbital flow is limited to the a narrow stream 
in the planet's neighbourhood, and the circumplanetary disc 
becomes more symmetric. This allows the planet to capture 
a larger amount of the mass, but this process is limited by 
the fact that the mass flux crossing the planet orbit decreases 
quickly with \Z\, thus ending the mass accumulation phase 
close to Z ~ —0.6. The fast mass accumulation builds up a 
strong pressure gradient in the planet's proximity, which is 
only supported by the gas inflow into the circumplanetary 
disc. After the mass accumulation phase ends gas starts to 
flow away from the circumplanetary disc, and short phase 
of outward migration can occur. 

The growth of the effective planet mass is strongly de- 
pendent on the assumed temperature profile in the circum- 
planetary disc, decreasing for larger circumplanetary disc 
aspect ratios. For h p > 0.5 we found no substantial mass 
accumulation, and gas is lost from the the planet's environ- 
ment during the entire slow migration limit. 

The planet mass plays an important role, even though 
the planet's orbital evolution is only weakly dependent on 
it in the fast migration limit. Since Z ~ Mp 2 ^ 3 , the planet 
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with bigger mass (or lower h p ) makes the transition to the 
slow migration limit sooner, and settle at larger radii than 
planets with lower mass (or larger h p ). 

The position where the planet stops rapid type III mi- 
gration thus depends on the planet's mass, but also on the 
disc structure. The mass of the co-orbital region defines the 
region of the disc where the rapid migration is allowed. To 
start and maintain type III migration the mass of the co- 
orbital region has to be comparable to the planet mass. This 
means that an increase of the disc mass gives a similar effect 
as a decrease of the planet's mass. 

The planet's orbital evolution is determined by the 
changes of the co-orbital mass deficit Ma, expressing the 
mass and the asymmetry of the co-orbital region. The mass 
is the product of the density profile and the volume of the 
co-orbital region. For discs with as > —2 it therefore de- 
creases with radius, and inward rapid migration always slows 
down. This is the stopping mechanism for most of the simu- 
lations presented in this paper. Type III migration can also 
be stopped at significant density jumps in the disc (such as 
a disc edge). In this case the asymmetry of the horseshoe 
region is lost due to the interaction with the low density re- 
gion and the planet rapidly makes the transition to the slow 
migration regime \Z\ <C 1. 

In all our simulations the eccentricity did not exceed 
0.04 and was damped during the rapid migration phase. 
We found the damping time-scale to be independent of 
the planet mass, but inversely proportional to the total 
disc mass. The measured damping time-scales are few times 
longer than the one predicted by the analytical formulation. 

The results in this paper show that type III migration 
can operate effectively, but does require the mass in the co- 
orbital region to be larger than the planet mass in order 
to get started, making it easier to start type III migration 
for lower mass planets, and for positions further out in the 
disk. The mass of the co-orbital region also sets a limit on 
how far the planet can migrate inward. The minimum orbital 
radius it can attain is again set by the position where the co- 
orbital mass exceeds twice the planet mass. It can however 
stop at larger radii than that, depending on the details of 
the density structure of the co-orbital region. For all cases 
the use of the dimensionless migration rate Z is useful to 
characterise the migration behaviour and the structure of 
the co-orbital region. Calculating Z can for example be used 
to check whether type I or type II migration will lead to type 
III migration. 

The simulations presented in this paper were set up to 
study inward migration in the type III regime. In a real 
disc-planet interaction, type III migration may also lead to 
outward migration, this depending on the details of the pre- 
vious migration pattern. This outward migration mode will 
be studied in detail in a subsequent paper. 
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